Electronic signature of the instantaneous asymmetry in the first coordination shell of 

liquid water 
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Interpretation of the X-ray spectra of water as evidence for its asymmetric structure has challenged 
the conventional symmetric nearly-tetrahedral model and initiated an intense debate about the order 
and symmetry of the hydrogen bond network in water. Here, we present new insights into the nature 
of local interactions in water obtained using a novel energy decomposition method. Our simulations 
reveal that while a water molecule forms, on average, two strong donor and two strong acceptor 
bonds, there is a significant asymmetry in the energy of these contacts. We demonstrate that 
this asymmetry is a result of small instantaneous distortions of hydrogen bonds, which appear as 
fluctuations on a timescale of hundreds of femtoseconds around the average symmetric structure. 
Furthermore, we show that the distinct features of the X-ray absorption spectra originate from 
molecules with high instantaneous asymmetry. Our findings have important implications as they 
help reconcile the symmetric and asymmetric views on the structure of water. 
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A detailed description of the structure and dynam- 
ics of the hydrogen bond network in water is essential 
to understand the unique properties of this ubiquitous 
and important liquid. It has long been accepted that 
the local structure of water at ambient conditions is 
tetrahedral [H-Q. Although the thermal motion causes 
distortions from the perfectly tetrahedral configuration, 
each molecule in the liquid is bonded, on average, to 
four nearest neighbours via two donor and two acceptor 
bonds [2j|. This traditional view is based on results from 
X-ray and neutron diffraction experiments 04^1, vibra- 
tional spectroscopy macroscopic thermodynam- 
ics data I2|. l3l llll as well as molecular dynamics simula- 
tions 0, mull • 

But this traditional picture has recently been ques- 
tioned based on data from the X-ray absorption, X-ray 
emission and X-ray Raman scattering experiments [HJ- 
Il9j . The results of these spectroscopic studies have been 
interpreted as evidence for strong distortions in the hy- 
drogen bond network with highly asymmetric distribu- 
tion of water molecules around a central molecule. It 
has been suggested that a large fraction of molecules 
form only two strong hydrogen bonds: one acceptor and 
one donor bond However, the "rings and chains" 

structure of liquid wat er Il6| . as well as the inhomoge- 
neous two-state model [17], [22| implied by such an inter- 
pretation, have been challenged on many fronts |3, M, |23| — 
31 1 and are a matter of an ongoing debate 0, SUSIU- 

In this work, we performed an unprecedented com- 
putational study of the energetics and symmetry of lo- 
cal interactions between water molecules in the liquid 
phase by utilizing a novel energy decomposition analysis 
based on absolutely localized molecular orbitals (ALMO 
EDA) [3a |. The decomposition of the interaction en- 
ergy into physically meaningful components provides a 
deeper insight into the nature and mechanisms of inter- 
molecular bonding than the traditional total-energy elec- 



tronic structure methods The success of ALMO 

EDA in investigating chemical bonding in molecular gas- 
phase complexes including small water clusters 45H47 
has inspired us to generalize this methodology and de- 
velop the first energy decomposition method for periodic 
condensed phase systems. 

ALMO EDA separates the total interaction energy of 
molecules (AEtot) into the interaction energy of the un- 
relaxed electron densities on the molecules (AEfrz) and 
the orbital relaxation energy. The latter can be further 
decomposed into an intramolecular relaxation associated 
with polarization of the electron clouds on molecules 
in the field of each other (AEpol), two-body donor- 
acceptor orbital interactions {AEdel) and a generally 
small higher-order (AEho ) relaxation term (see Ref. [38| 
for a detailed description of the ALMO EDA terms). 

AEtot = AEfrz + AEpol + AEdel + AEho 

Mol Mol 

AEdel = ]T AE C = &E d ^a (1) 
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Two-body components AEd^a are the main focus of 
this work and arise due to derealization of electrons 
(or charge transfer) from the occupied orbitals of donor 
molecule D to the virtual orbitals of acceptor A. Each 
of these terms provide an accurate measure of the per- 
turbation of orbitals localized on a molecule by donor 
or acceptor interactions with a single neighbour in the 
many-body system. The donor-acceptor energies are ob- 
tained self-consistently and include cooperativity effects, 
which are essential for a correct description of the hy- 
drogen bond networks (U [48|. Furthermore, the two- 
body terms are natural local descriptors of intermolecular 
bonding, which allowed us to analyze the molecular net- 
work in liquid water without introducing any arbitrary 
definitions of a hydrogen bond [49| . 

The unique insight obtained by applying ALMO EDA 
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to liquid water revealed a significant asymmetry in the 
strength of the local donor-acceptor contacts. Ab initio 
molecular dynamics (AIMD) simulations performed with 
the recently developed second-generation Car-Parrinello 
approach [50| enabled us to characterize geometric ori- 
gins of the asymmetry, its dynamical behaviour, as well 
as the mechanism of its relaxation. Furthermore, to ad- 
dress the controversial question of whether it is correct 
to interpret the X-ray spectra of water in terms of asym- 
metric structures we performed extensive calculations of 
its X-ray absorption (XA) spectrum and compared the 
spectral characteristics of water molecules with different 
degree of asymmetry. 

RESULTS 

Electronic asymmetry and its origins. ALMO 
EDA was performed for the decorrclated configurations 
collected from an AIMD simulation at constant temper- 
ature (300 K) and density (0.9966 g/cm 3 , see Methods 
for computational details). The electron derealization 
energy per molecule AEc defined in Equation ([1]) was 
analysed by considering each water molecule as a donor 
and as an acceptor: 

Mol Mol 

AEc = ^ AEc^n = ^ AEpf^c, 

N=l iV=l 

where C stands for the central molecule and N for its 
neighbours. It is important to emphasize that terms 
donor and acceptor are used throughout this paper to 
describe the role of a molecule in the transfer of the elec- 
tron density. This is opposite to the labeling used for a 
donor and an acceptor of hydrogen in a hydrogen bond. 

Figure Q] shows contributions of the five strongest 
donor-acceptor interactions to the average derealization 
energy of a molecule (AEc) (brackets (..) denote averag- 
ing over all central molecules and AIMD snapshots). It 
demonstrates that electron derealization is dominated 
by two strong interactions, which together are respon- 
sible for ^93% of the derealization energy of a single 
molecule. The third and the fourth strongest donor (ac- 
ceptor) interactions contribute ^5% and correspond to 
back-donation of electrons to (from) the remaining two 
first-shell neighbours (i.e. there is non-negligible der- 
ealization from a typical acceptor to a typical donor). 
The remaining ^2% correspond to the derealization en- 
ergy to (from) the second and more distant coordination 
shells. 

Comparison of the strengths of the first and second 
strongest donor-acceptor interactions (~25 kJ/mol and 
~12 kJ/mol, respectively) with that in water dimcr 
(~9 kJ/mol) suggests that each water molecule can be 
considered to form, on average, two donor and two ac- 
ceptor bonds. Substantial difference in the strengths of 
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FIG. 1. Five strongest interactions. Average con- 
tributions of the strongest acceptor (AEn^c) and donor 
(AEc^n) interactions. The rightmost column shows that 
higher-order derealization (AEho) does not contribute to 
the overall binding significantly. 

the first and second strongest interactions implies that a 
large fraction of water molecules experience a significant 
asymmetry in their local environment. To characterize 
the asymmetry of the two strongest donor contacts of a 
molecule we introduced a dimensionlcss asymmetry pa- 
rameter 

AE c ^N 2 " d 
D = ~AE 

and an equivalent parameter for the acceptor interactions 
Ta- The asymmetry parameter is zero if the two con- 
tacts are equally strong and equals to one if the second 
contact does not exist. The probability distribution of 
molecules according to their T-parameters is shown in 
Figure [2] together with the lines separating the molecules 
into four groups of equal sizes with different asymme- 
try. The shape of the distribution demonstrates that 
most molecules form highly asymmetric donor or accep- 
tor contacts at any instance of time. For example, the 
line at T ps 0.5 indicates that for ~75% of molecules ei- 
ther Ta or T u is more than 0.5, which means that the 
strongest donor or acceptor contact is at least two times 
stronger than the second strongest for these molecules. 
Furthermore, the peak in the region of high T in Figure [2] 
indicates the presence of molecules with significantly dis- 
torted hydrogen bonds. 

To understand the origins of the asymmetry we com- 
pared the geometry of donor-acceptor pairs involved in 
the first and second strongest interactions. We found 
that the strength of the interaction is greatly affected 
by the intermolecular distance R = d(On — Oa) and 
the hydrogen bond angle (3 = /LOdOaH while the 
other geometric parameters have only minor influence 
on AEd^a- The strongly overlapping distributions in 
Figure [3] suggest that some second strongest interactions 
have the same energetic and geometric characteristics as 
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FIG. 2. The normalized probability density function of 
the asymmetry parameters Ta and Td. The probability 
of finding a molecule in a bin can be found by dividing the 
corresponding density value by the number of bins (i.e. 400). 
The dashed black lines at Ta | , | , | partition all molecules 
into four groups of equal sizes. 



the strongest contacts. This implies that the observed 
electronic asymmetry cannot be attributed to the pres- 
ence of two distinct types of hydrogen bonds - weak and 
strong. It is rather a result of continuous deformations 
of a typical bond. Another important conclusion that 
can be made from the distributions in Figure |3] is that 
relatively small variations of the intermolecular distance 
(R ~ 0.2 A) and hydrogen bond angle (/3 ~ 5 - 10°) lead 
to the remarkable changes in the strength and electronic 
structure of hydrogen bonds. Analysis of the structure of 
the molecular chains defined by the first strongest bonds 
(i.e. one donor and one acceptor for each molecule) shows 
that their directions are random, without any long-range 
order (i.e. rings, spirals or zig-zags) on the length scale 
of the simulation box (~ 15 A). 

We did not attempt to quantify the concentration of 
single-donor and single-acceptor molecules in this work 
because the slow decay of the distribution tails (Fig- 
ure [3]) implies that defining such configurations using a 
distance, angle or energy cutoff is an unavoidably arbi- 
trary procedure. A quantitative analysis of the network, 
which was performed for the same molecular configura- 
tions in Ref. Il5l shows that according to the commonly 
used geometric definitions of hydrogen bonds [l6|, |j3l], HH 
the structure of water is distorted tetrahedral with only 
a small fraction of broken bonds. The results presented 
in Figure [3] indicate that these geometric criteria cannot 
fully characterize the dramatic effect of distortions on the 
local electronic structure and donor-acceptor interactions 
of water molecules. 

It is important to note that some second strongest in- 
teractions are weakened by distortions to such an extent 
that back-donation to (from) a nearby donor (acceptor) 
becomes the second strongest interaction. Such config- 
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FIG. 3. Energetic and geometric characteristics of the 
instantaneous asymmetry. The probability distribution 
of the (a) strength, (b) intermolecular distance R and (c) 
hydrogen bond angle j3 for the first (red) and second (blue) 
strongest donor interactions C — > N. Filled areas show the 
contribution of configurations, for which back-donation to a 
nearby donor is stronger than donation to the second acceptor 
(cf. text). Distributions for acceptor interactions, N — > C, 
are almost identical and not shown. 



urations can be clearly distinguished by the large-angle 
peak in Figure (the region of the 10-fold magnifica- 
tion). They account for ~ 6 — 7% of all configurations 
and are responsible for the low-energy peak in the distri- 
bution of AEd^a (filled blue areas in Figure [3]) and for 
the high-T peak in Figured! 

Relaxation of the instantaneous asymmetry. 

The overlapping distributions in Figure[3]suggest that de- 
spite the difference in the strength of the donor-acceptor 
contacts their nature is similar and in the process of ther- 
mal motion the strongest interacting pair can become the 
second strongest and vice versa. To estimate the time 
scale of this process, we examined how the average energy 
of the first two strongest interactions fluctuates in time 
(see Methods for definitions of time-dependent averages). 
Figure @Ji shows that the strength of these interactions 
oscillates rapidly and after ^80 fs from an arbitrarily cho- 
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sen time origin, the first strongest interaction becomes 
slightly weaker than the second strongest (note that first 
and second refer to their order at t = 0). The amplitude 
of the oscillations decreases and the strengths of both in- 
teractions approach the average value of ^20 kJ/mol on 
the time scale of hundreds of femtoseconds. The decay of 
the oscillations indicates fast decorrelation of the time- 
separated instantaneous values because of the strong cou- 
pling of a selected pair of molecules with its surround- 
ings. In other words, although the energy of a particular 
hydrogen bond fluctuates around its average value in- 
definitely (i.e. with a never-decreasing amplitude), this 
bond has approximately equal chances of becoming weak 
or strong after a certain period of time independently 
of its strength at t = 0. This effect is due to the noise 
introduced by the environment and can be observed in 
ultrafast infrared spectroscopy experiments Q • 

It is important to emphasize that the time averages 
shown in Figure |4] are physically meaningful and can 
be calculated accurately only for the time intervals that 
are shorter than the average lifetime of a hydrogen bond 
thb ~ 5 ps [TBI, 51 1 (see Methods). The small residual 
asymmetry that is still present after 500 fs (Figure 2£i) is 
an indication of the slow non-exponential relaxation be- 
haviour that characterises the kinetics of many processes 
in liquid water [5~0 ]. 

In addition to the instantaneous values of AEd^a oi 
time t, the dashed lines in Figure^ show the correspond- 
ing averages over time t (see Methods) . The latter shows 
that any neighbour-induced asymmetry in the electronic 
structure of a water molecule can be observed only with 
an experimental probe with a time-resolution of tens of 
femtoseconds or less. On longer time scales the asymme- 
try is destroyed by the thermal motion of molecules and 
only the average symmetric structures can be observed 
in experiments with low temporal resolution. 

An examination of time dependence of all two-body 
and some three-body geometric parameters that charac- 
terize the relative motion of molecules reveals the mech- 
anism of the relaxation. Similar shapes of the curves in 
Figures and St shows that the relaxation of the asym- 
metry is primarily caused by low-frequency vibrations of 
the molecules relative to each other. The minor differ- 
ences in the behaviour of the curves, in particular at 80 fs, 
indicate that the relaxation of the asymmetry is also in- 
fluenced by some other degrees of freedom. The temporal 
changes in the hydrogen bond angles towards the aver- 
age value (Figure Ht) show that librations of molecules 
play this minor role in the relaxation process. The ab- 
sence of the correlation between the energy curves and 
the remaining geometric parameters indicates that the 
other types of motion are not important for the relax- 
ation mechanism. 

The kinetics and mechanism of the asymmetry relax- 
ation presented here are supported by data from ultrafast 
infrared spectroscopy, which can directly observe inter- 
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FIG. 4. Relaxation of the instantaneous asymme- 
try. Time dependence of the (a) average strength, (b) in- 
termolecular distance R and (c) hydrogen bond angle j3 for 
the first (red) and second (blue) strongest donor interactions 
C — i> N . Solid lines show the instantaneous values while 
the dashed lines correspond to the time-average values (see 
Methods). Time-dependent characteristics of acceptor inter- 
actions, N —¥ C, are almost identical and not shown. 



molecular oscillations with a period of 170 fs . They are 
also in agreement with the theoretical work of Artacho 
et al., who have used the Mulliken bond order parameter 
to characterise the connectivity and dynamical processes 
in the hydrogen bond network of liquid water [29j | . 

X-ray absorption signatures of asymmetric 
structures. The time behaviour described above im- 
plies that the instantaneous asymmetry can, in princi- 
ple, be detected by X-ray spectroscopy, which has tem- 
poral resolution of several femtoseconds and is highly 
sensitive toperturbations in the electronic structure of 



molecules [16|, l21| . To identify possible relationships be- 
tween the spectroscopic features and asymmetry, we cal- 
culated the XA spectrum of liquid water at the oxygen K- 
edge using the half-core-hole transition potential formal- 
ism 53] within all-electron density functional theory (3 



(see Methods for computational details). Although the 
employed computational approach overestimates intensi- 
ties in the post-edge part of the spectrum and underesti- 
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mates the pre-edge peak and overall spectral width [55| 
it provides an accurate description of the core-level exci- 
tation processes and semiquantitatively reproduces the 
main features of the experimentally measured spectra 
(Figure [Sh). 

The localized nature of the Is core orbitals allows to 
disentangle spectral contributions from molecules with 
different asymmetry. To this end, all molecules were sep- 
arated into four groups according to the asymmetry of 
their donor and acceptor environments as shown in Fig- 
ure [2] Choosing boundaries at T w | , | , | distributes 
all molecules into four groups of approximately equal 
sizes (i.e. 25 ± 2%). Figure [5jj shows four XA spec- 
tra obtained by averaging the individual contributions 
of molecules in each group. It reveals that molecules 
in the symmetric environment exhibit strong post-edge 
peaks while molecules with high asymmetry of their en- 
vironment are characterized by the amplified absorption 
in the pre-edge region. Furthermore, the relationship be- 
tween the asymmetry and absorption intensity is nonuni- 
form: the pre-edge peak is dramatically increased in the 
spectrum for the 25% of molecules in the most asym- 
metric group, for which the first strongest interaction is 
more than six times stronger than the second. As a con- 
sequence, the pre-edge feature of the calculated XA is 
dominated by the contribution of molecules in the highly 
asymmetric environments (Figure [5J;) . 

The pronounced pre-edge peak in the experimentally 
measured XA spectrum of liquid water has been inter- 
preted as evidence for its "rings and chains" structure, 
where ^80% of molecules have two broken hydrogen 
bonds 1(3, 2(j. Our results suggest that this feature of 
the XA spectrum can be explained by the presence of 
a smaller fraction of water molecules with high instan- 
taneous asymmetry. Although the employed XA mod- 
eling methodology does not allow us to estimate pre- 
cisely the size of this fraction our result is consistent with 
that of recent theoretical studies at an even higher level 
of theory, which have demonstrated that the main fea- 
tures of the experimental XA spectra can be reproduced 
in simulations based on conventional nearly-tetrahcdral 
models (3ll . 56 1. Our work complements the previous re- 



sults by revealing an interesting and important connec- 
tion between relatively small geometric perturbations in 
the hydrogen-bond network, the large asymmetry in the 
electronic ground state and the XA spectral signatures 
of the core-excitation processes. 

Conclusions. We have applied a combination 
of ALMO EDA, XA spectra modelling and second- 
generation Car-Parrinello molecular dynamics, to study 
perturbation induced in the electronic structure of water 
molecules by local donor-acceptor interactions with their 
neighbours in the liquid phase. Our main conclusions 
are as follows. First, the strength of donor-acceptor in- 
teractions suggest that each molecule in liquid water at 
ambient conditions forms, on average, two donor and two 
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FIG. 5. X-ray absorption (XA) spectra of liquid wa- 
ter, (a) Comparison of the calculated and experimental 
XA spectra, (b) Calculated XA spectra of the four groups 
of molecules separated according to the asymmetry of their 
donor (Td) and acceptor (Ta) environments as shown in the 
inset, (c) Contributions of the four groups into the total XA 
spectrum. The color coding is shown in the inset above. 



acceptor bonds. Second, even small thermal distortions 
in the tetrahedral hydrogen bond network induce a sig- 
nificant asymmetry in the strength of the contacts caus- 
ing one of the two donor (acceptor) interactions to be- 
come, at any instance of time, substantially stronger than 
the other. Thus, the instantaneous structure of water is 
strongly asymmetric only according to the electronic cri- 
teria, not the geometric one. Third, intermolecular vibra- 
tions and librations of OH groups of hydrogen bonds re- 
sult in the relaxation of the instantaneous asymmetry on 
the time scale of hundreds of femtoseconds. Finally, the 
pronounced pre-edge peak observed in the XA spectra of 
liquid water can be attributed to molecules in asymmet- 
ric environments created by instantaneous distortions in 
the fluctuating but on average symmetric hydrogen bond 
network. 

Thus, our work helps reconcile the two existing and 
seemingly opposite models of liquid water - the tradi- 
tional symmetric and the recently proposed asymmetric 
- and represents an important step towards a better un- 
derstanding of the electronic structure of the hydrogen 
bond network of one of the most important liquids on 
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Earth. 

The broader significance of this work lies in the de- 
velopment of an energy decomposition method for con- 
densed phase systems, which can provide deep insight 
into the electronic structure of a wide range of molecular 
systems and, thus, offers new opportunities for studying 
their complex behaviour. 



METHODS 

Ab initio molecular dynamics. Ah initio MD simu- 
lations with classical nuclei based on the Pcrdcw-Burkc- 
Ernzerhof (PBE) exchange-correlation (XC) functional 
were performed at constant temperature (300 K) and 
density (0.9966 g/cm 3 ). Long 70 ps AIMD trajectories 
for systems containing 128 water molecules were obtained 
by taking advantage of the computational efficiency of 
the second-generation Car-Parrinello method as de- 
scribed in Ref. 15. Our previous work fl5l ] has shown 
that PBE provides a realistic description of many im- 
portant structural and dynamical characteristics of liq- 
uid water including the radial distribution functions, self- 
diffusion and viscosity coefficients, vibrational spectrum 
and hydrogen-bond lifetime. However, the PBE water 
with the classical description of the nuclei is somewhat 
overstructured suggesting that the degree of the network 
distortion, fraction of broken bonds and, thus, the asym- 
metry in our simulations may be underestimated. There- 
fore, we verified that the main conclusions presented in 
this work are also valid for snapshots generated with sim- 
ulations with quantum hydrogen nuclei and based on a 
water model, which rectifies the main shortcomings of 
the PBE functional (sec Supplementary Information). 

Energy decomposition analysis. ALMO EDA [38[ 
for periodic systems was implemented [57j in the CP2K 
package, which relies on the mixed Gaussian and plane 
wave (GPW) representation of the electronic degrees of 
freedom [H^]. The GPW approach makes CP2K uniquely 
suited for ALMO EDA since the localized atom-centered 
Gaussian basis sets are required for the construction of 
absolutely localized molecular orbitals, whereas the use 
of plane waves ensures computational efficiency of the 
large-scale calculations performed in this work. 

Molecular orbitals in all ALMO EDA calculations were 
represented by a triple-^ Gaussian basis set with two sets 
of polarization functions (TZV2P) optimized specifically 
for molecular systems [5^. A high density cutoff of 1000 
Ry was used to describe the electron density. The XC en- 
ergy was approximated with the BLYP functional. The 
Brillouin zone was sampled at the L-point and separable 
norm-conserving pseudopotentials were used to describe 
the interactions between the valence electrons and the 
ionic cores |(3(| ■ We verified that performing ALMO EDA 
with the HSE06 screened hybrid XC functional, which 
provides better description of band gaps and electron 



derealization effects than BLYP, does not change the 
main conclusions of the work (see Supplementary Infor- 
mation). 

Averaging procedures. The average donor-acceptor 
energies shown in Figure Q] and distributions in Figures [5] 
and [3] were obtained by performing ALMO EDA for 701 
decorrelated AIMD snapshots separated by 100 fs (i.e. 
89,728 molecular configurations). 

The time-dependent averages in Figure |4jt were com- 
puted using the ALMO EDA terms for 3,501 snapshots 
separated by 20 fs (448,128 local configurations). The 
instantaneous values at time t (solid lines in Figure [4^) 
were calculated by averaging over time origins r sepa- 
rated by 100 fs and over all surviving triples: 

l t 1 M(r,t) 
(AE c ^ N (t)) = -^——r ]T AE C ^ N (T + t), 

T=l ' ' C=l 

where M(t, t) is the number of triples that survived from 
time r to r+t. A triple is considered to survive a specified 
time interval if the central molecule has the same two 
strongest-interacting neighbours in all snapshots in this 
interval. 

The average values over time t (dashed lines in Fig- 
ure^) were calculated by averaging over time origins r, 
all snapshots lying in the time interval from r to r + 1 
and over all surviving triples: 

(AE c ^ N (t))* = 

T t M{r,t) 

= -Y—Y—!-^ y ae c ^ n (t+k) 

T t— 1 1 + 1 M(t, t) 

r =l K =0 v ' 1 C=l 

Time-averages (A(t))* are related to the instantaneous 
values (A(t)) by the following equation: 

(A(t)r^- t j\A(r))dT, 

where equality holds if all triples survive over time t. 

The time averages defined here are physically mean- 
ingful and can be accurately calculated only for the time 
intervals that are shorter than the average lifetime of 
a triple or, equivalently, the average lifetime of a hy- 
drogen bond, t H b ~ 5 ps 0, 51 1. Nevertheless, some 
triples break apart even after shorter time affecting the 
time averages. For example, only ^80% of all triples sur- 
vive over the time interval shown in Figure |4j We found 
that stronger interacting triples have longer lifetimes and, 
therefore, are overrepresented at long time. However, the 
error introduced by this bias is only ~1.5 kJ/mol at 0.5 ps 
and does not affect the main conclusions. This error is 
the main reason that the long-time averages in Figure [4jt 
converge to 20 kJ/mol instead of the expected value of 
18.5 kJ/mol. 

X-ray absorption spectrum. The XA spectra of 
water at the oxygen K-edge were obtained using the all- 
electron Gaussian augmented plane wave formalism im- 
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plemented in the CP2K package [54], 155J. The calcula- 
tions were performed with the BLYP XC functional and 
half-core-hole transition potential. Large basis sets (6- 
31 1G** for hydrogen and cc-pVQZ for oxygen atoms) 
were used to provide an adequate representation of the 
unoccupied molecular orbitals in the vicinity of absorbing 
atoms [55|] . A density cutoff of 320 Ry was used to de- 
scribe the soft part of the electron density. The onset en- 
ergies of the absorption spectra were aligned with the first 
ASCF excitation energy obtained for each oxygen atom 
from a separate calculation with the same setup. The 
spectral intensities were calculated as transition moment 
integrals in the velocity form. To mimic the experimen- 
tal broadening, the discrete lines were convoluted with 
Gaussian functions of 0.2 cV width at half- maximum. 
The final spectrum was obtained by averaging the con- 
voluted spectra of 9,024 oxygen atoms from 141 AIMD 
snapshots separated by 500 fs (i.e. 64 oxygen atoms were 
randomly selected in each snapshot). 
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